%% 基于dogleg方法的信赖域算法Trust_region_search 
% 初始化
eps=1e-8;% 梯度阈值判断，当计算点处梯度小于eps时结束算法
lb=[-100,-100];% 变量最小值
ub=[100,100];% 变量最大值
n=length(lb);% 变量个数
rand1=rand();
x0=rand1*ub+(1-rand1)*lb;% 迭代初始点
maxiter=1000;% 最大迭代次数
eta1=0.25;% rho值判断：小于eta1时说明需减小信赖域半径
eta2=0.75;% rho值判断：大于eta2时说明可以增大信赖域半径
delta=50;
delta_max=100;% 最大信赖域半径
eta=0;% 初始化eta
x=x0;% 初始化x点位置
i=0;% 初始化迭代次数
while  i < maxiter
    g = gradient(x);% 求x点梯度
    if norm(g)<= eps
        break;
    end
    h = hessian(x);% 海瑟矩阵或其近似矩阵
    p = subproblem(g,h,delta);% 利用dogleg方法求得搜索向量p
    rho =-(f(x)-f(x+p))/(g'*p+0.5*p'*h*p);% 计算rho值
    if rho < eta1 % 该部分简单判断rho值并对信赖域半径delta判断
        delta = eta1*delta;
    elseif rho>eta2 && abs(norm(p)-delta)<1e-8
        delta = min(2*delta,delta_max);
    end
    if rho>eta
        x = x+p;
    end
    i = i+1;
end
y = x;% y值为最终求得点
%% Dogleg方法求解子问题
function p = subproblem(g,h,delta)
p_B = -inv(h)*g;% 在牛顿法搜索方向上极小值点位置向量
p_U = -g'*g/(g'*h*g)*g;% 在梯度下降方向上极小值点位置向量
if norm(p_B)<= delta %当牛顿法搜索方向向量小于信赖域半径时，返回该向量
    p = p_B;
elseif norm(p_U) >= delta % 当梯度下降方向向量大于信赖域半径时，取该方向，
    % 步长为信赖域半径
    p = delta/norm(p_U)*p_U;
else % 不满足以上情况时，考虑两方向综合
    p_BU = p_B-p_U;
    tau1 = (-p_U'*p_BU+sqrt((p_U'*p_BU)^2-p_BU'*p_BU*(p_U'*p_U-delta^2)))/(p_BU'*p_BU);
    if tau1>=1 && tau1<=2
        tau = tau1;
    else
        tau = (-p_U'*p_BU+sqrt((p_U'*p_BU)^2-p_BU'*p_BU*(p_U'*p_U-delta^2)))/(p_BU'*p_BU);
    end
    p = p_U+tau*p_BU;
end
end
% 目标函数
function f = f(x)
f = 100*(x(1)^2 - x(2))^2 + (x(1) - 1)^2;
end
% 目标函数梯度
function g = gradient(x)
dg_dx1 = 2*(x(1)-1)+400*x(1)*(x(1)^2-x(2));
dg_dx2 = 200*(x(2)-x(1)^2);
g = [dg_dx1;dg_dx2];
end
% 目标函数海瑟矩阵或其近似矩阵
function h = hessian(x)
h_11 = 1200*x(1)^2-400*x(2)+2;
h_12 = -400*x(1);
h_21 = h_12;
h_22 = 200;
h = [h_11,h_12;h_21,h_22];
end